#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <string.h> // memcpy
#include <fstream>
/*
  use M_PI const: 3.14159 
*/

#define A 0
#define B 2
#define SIGMA 1

double psi(double x){ // функция, задающая начальное распределение электронной плотности
  //return ((x >= 0) && (x <= 0.5)) ? cos(2*M_PI*x)*sin(2*M_PI*x) : 0;
  return ((x >= 0) && (x <= 0.5)) ? sin(2*M_PI*x) : 0;
}

double fi_a(double time){ // левое граничное условие
  return 0;
}
double fi_b(double time){ // правое граничное условие
  return 0;
}
double f(double currU, double time, double x){ // функция, задающая внешнее воздействие
  return (((x >= 0.58) && (x <= 0.6))||((x >= 0.78) && (x <= 0.8))) ? 5 : 0; //2 : 0;
  //return ((x >= 0) && (x <= 2)) ? 2 : 0;
}
double getDeltaT(double deltaX){
  return (deltaX * deltaX) / (2 * SIGMA);
}
int getStepsT(double deltaT, double time){
  return (int)((time / deltaT) + 1);
}
int getStepsX(double deltaX){
  return (int)(((B - A) / deltaX) + 1);
}

double **calculation(
                      double deltaX, // шаг по X
		     double time, // кол-во моделируемого времени
		      double deltaTShow // шаг во времени, с которым запоминаем
		    ){
  // считаем все, запоминаем нужное
  
  // шаг по времени определяется из ограничений
  double deltaT = getDeltaT(deltaX);
  
  double **data;
  int stepsT = getStepsT(deltaT, time);
  int stepsX = getStepsX(deltaX);
  data = new double*[(int)(time/deltaTShow) + 1]; // ряды в рассматриваемые моменты времени
  for(int i = 0; i < ((int)(time/deltaTShow) + 1); i++){
    data[i] = new double [stepsX]; // массивы значений в момент времени
  }
  
  double **temp;
  temp = new double*[2]; // на смену 2 временных ряда
  temp[0] = new double [stepsX];
  temp[1] = new double [stepsX];
  
  // считаем:
  data[0][0] = fi_a(0); // левое граничное значение плотности в нулевой момент времени
  data[0][stepsX - 1] = fi_b(0); // правое граничное значение плотности в нулевой момент времени
  for(int i = 1; i < (stepsX - 1); i++){ // плотность распределения в нулевой момент времени
    data[0][i] = psi(deltaX * i);
  }
  memcpy(temp[0], data[0], sizeof(data[0]) * stepsX);
  int iShow = 1;
  for(int i = 1; i < stepsT; i++){ // строим электронные плотности во все моменты времени
    //data[i][0] = fi_a(i); // левое граничное значение плотности в i-й момент времени
    temp[i%2][0] = fi_a(i);
    //data[i][stepsX - 1] = fi_b(i); // правое граничное значение плотности в i-и момент времени
    temp[i%2][stepsX - 1] = fi_b(i);
    for(int j = 1; j < (stepsX - 1); j++){
      temp[i%2][j] = 
         temp[(i-1)%2][j] // значение плотности в данной точки в прошлый момент времени
       + ((SIGMA * deltaT)/(deltaX * deltaX))
       * (
	        temp[(i-1)%2][j-1] // левое соседнее значение плотности в прошлый момент
	  - 2 * temp[(i-1)%2][j] // значение плотности в данной точке в прошлый момент
	  +     temp[(i-1)%2][j+1] // правое соседнее значение плотности в прошлый момент
         )
       + f( // внешнее воздействие, зависит от:
	    temp[(i-1)%2][j] // значение плотности в данной точке в прошлый момент времени
	   ,(i-1) * deltaT // значение времени на прошлом шаге
	   ,j * deltaX // координата данной точки
          ) * deltaT // шаг временной сетки
      ;
      //std::cout << i*deltaT << "\n";
     // std::cout << () << "\n";
      
      if( i*deltaT > deltaTShow*iShow){ // как только, так сразу запоминаем
        
	memcpy(data[iShow], temp[i%2], sizeof(temp[0]) * stepsX);
	iShow++;
	//std::cout << iShow << "\n";
      }
    }
  }
  return data;
}

int main(int argc, char **argv){ 
  double deltaX = atof(argv[1]); // шаг сетки по х // в статье h
  double time = atof(argv[2]); // время моделирования
  double deltaTShow = atof(argv[3]); // шаг по времени в результате

  
  double **data = calculation(deltaX, time, deltaTShow);
  double deltaT = getDeltaT(deltaX);
  int stepsT = getStepsT(deltaT, time);
  int stepsX = getStepsX(deltaX);
  /*
  int deltaTOutCoeff = (deltaT < deltaTShow) ? ((int)(deltaTShow / deltaT)) : 1;
  int stepsTOut = (int)(stepsT / deltaTOutCoeff);
  double **dataOut;
  dataOut = new double*[stepsTOut]; // ряды в моменты времени
  for(int i = 0; i < stepsTOut; i++){
    dataOut[i] = data[i*deltaTOutCoeff]; // массивы значений в момент времени
  }
  */
  /*std::cout << "Рабочий отрезок [" << A << ", " << B << "]\n";
  std::cout << "SIGMA: " << SIGMA << "\n";
  std::cout << "Шаг сетки по времени: " << deltaT << "\n";
  std::cout << "Шаг сетки по X: " << deltaX << "\n";
  std::cout << "Шагов по времени: " << stepsT << "\n";
  std::cout << "Шагов по X: " << stepsX << "\n";*/
  std::cout << "{\"u\":[";
  for(int i = 0; i < (int)(time/deltaTShow); i++){
    std::cout << "[";
    for(int j = 0; j < stepsX; j++){
      std::cout << data[i][j];
      if(j < (stepsX - 1)){
        std::cout << ",";
      }
    }
    std::cout << "]";
    if(i < ((int)(time/deltaTShow) - 1)){
      std::cout << ",";
    }
  }
  std::cout << "],\"f\":[";
    for(int i = 0; i < (int)(time/deltaTShow); i++){
    std::cout << "[";
    for(int j = 0; j < stepsX; j++){
      std::cout << f(data[i][j], i * deltaTShow, j * deltaX);
      if(j < (stepsX - 1)){
        std::cout << ",";
      }
    }
    std::cout << "]";
    if(i < ((int)(time/deltaTShow) - 1)){
      std::cout << ",";
    }
  }
  std::cout << "]}";
  return 0;
}